/*
  Copyright (C) 2011, Kyungjoo Kim
  
  This file is part of UHM (Unassembled HyperMatrix) sparse direct solver.
  
  All rights reserved.

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public
  License as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.
  
  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.
  
  You should have received a copy of the GNU Lesser General Public
  License along with this library; if not, write to the Free Software
  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  Also add information on how to contact you by electronic and paper mail.
  
  
  Kyungjoo Kim
  iamkyungjoo@gmail.com
*/
#include "uhm/common.hxx"
#include "uhm/const.hxx"

#include "uhm/matrix/uhm/matrix.hxx"
#include "uhm/matrix/uhm/fla.hxx"

namespace uhm {

  static inline int lu_nopiv_flat( int fs, int ss, 
				   linal::Flat_ ATL, linal::Flat_ ATR,
				   linal::Flat_ ABL, linal::Flat_ ABR );

  static inline int lu_nopiv_hier( int fs, int ss,
				   linal::Hier_ ATL, linal::Hier_ ATR,
				   linal::Hier_ ABL, linal::Hier_ ABR );
  
  
  void Matrix_FLA_::lu_nopiv() {

#ifdef UHM_HIER_MATRIX_ENABLE
    // ----------------------------------------------------------
    // ** Hier-Matrix
    // ----------------------------------------------------------
    lu_nopiv_hier( this->fs, this->ss, 
		   this->hier.ATL, this->hier.ATR,
		   this->hier.ABL, this->hier.ABR );
#else
    // ----------------------------------------------------------
    // ** Flat-Matrix
    // ----------------------------------------------------------
    lu_nopiv_flat( this->fs, this->ss, 
		   this->flat.ATL, this->flat.ATR,
		   this->flat.ABL, this->flat.ABR );
#endif
  }
  static inline int lu_nopiv_flat( int fs, int ss, 
				   linal::Flat_ ATL, linal::Flat_ ATR,
				   linal::Flat_ ABL, linal::Flat_ ABR ) {
    if (fs) 
      FLA_LU_nopiv(~ATL);
  
    if (fs && ss) {
    
#pragma omp task firstprivate(ATL, ATR)
      FLA_Trsm( FLA_LEFT, FLA_LOWER_TRIANGULAR,
		FLA_NO_TRANSPOSE, FLA_UNIT_DIAG,
		FLA_ONE, ~ATL, ~ATR );

#pragma omp task firstprivate(ATL, ABL)
      FLA_Trsm( FLA_RIGHT, FLA_UPPER_TRIANGULAR,
		FLA_NO_TRANSPOSE, FLA_NONUNIT_DIAG,
		FLA_ONE, ~ATL, ~ABL );

#pragma omp taskwait

      FLA_Gemm( FLA_NO_TRANSPOSE, FLA_NO_TRANSPOSE,
		FLA_MINUS_ONE, ~ABL, ~ATR, FLA_ONE, ~ABR );
    }
    return true;
  }

  static inline int lu_nopiv_hier( int fs, int ss, 
				   linal::Hier_ ATL, linal::Hier_ ATR,
				   linal::Hier_ ABL, linal::Hier_ ABR ) {

    if (fs) {
      linal::lu_nopiv(ATL);
 
      if (ss) {

#pragma omp task firstprivate(ATL, ATR)
        linal::trsm( FLA_LEFT, FLA_LOWER_TRIANGULAR,
                            FLA_NO_TRANSPOSE, FLA_UNIT_DIAG,
                            FLA_ONE, ATL, ATR );
        
#pragma omp task firstprivate(ATL, ABL)
        linal::trsm( FLA_RIGHT, FLA_UPPER_TRIANGULAR,
                            FLA_NO_TRANSPOSE, FLA_NONUNIT_DIAG,
                            FLA_ONE, ATL, ABL );
        
#pragma omp taskwait
        
        linal::gemm( FLA_NO_TRANSPOSE, FLA_NO_TRANSPOSE,
                            FLA_MINUS_ONE, ABL, ATR, FLA_ONE, ABR );

      }
    }
    return true;
  }
 
}
